{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "K-means clustering\n",
    "===========\n",
    "\n",
    "The k-means clustering algorithm is an unsupervised learning method broadly used in cluster analysis. the algorithm is based on two main steps:\n",
    "\n",
    "1. **Assignment**: each observation in the input space is assigned to the Best Matching Unit (BMU) based on a distance measure.\n",
    "2. **Update**: for each cluster estimate the mean and assign it to the centroids.\n",
    "\n",
    "In general the distance measure used is th [Euclidean distance](https://en.wikipedia.org/wiki/Euclidean_distance). A more efficient version is the [cosine distance](https://en.wikipedia.org/wiki/Cosine_similarity) that is based on the inner product. The goal in k-means clustering is to minimize the within-cluster sum of squares (reconstruction error):\n",
    "\n",
    "$$\\underset{S}{\\mathrm{argmin}} \\sum_{i=1}^{k} \\sum_{\\bf{x} \\in S_{i}} \\lVert \\bf{x} - \\bf{\\mu}_{i} \\rVert^{2}$$\n",
    "\n",
    "where $\\bf{x}$ are the input vectors, $S = \\{ S_{1},...,S_{k} \\}$ are the number of clusters, and $\\mu$ is the mean of the vectors belonging to each cluster."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Online k-means algorithm in Tensorflow\n",
    "--------------------------------------------\n",
    "\n",
    "Here I show you how to create an online version of the k-means algorithm. By *online* I mean that a batch of examples can be passed every time the training operation is called. This is particularly useful if the input space is particularly big. First of all we declare the input size and the number of clusters we want to use:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "input_size = 6\n",
    "number_centroids = 2 #the number of centroids\n",
    "batch_size = 3 #The number of input vector to pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can write the real algortihm. The centroids and the input are assigned through two placeholders. The distance used is the **Euclidean distance**. Remember that we want to minimize the within-cluster sum of squares (reconstruction error) between the input vectors and the centroids. Since the problem can be defined in term of a loss function, we can use our dear optimizers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import tensorflow as tf\n",
    "\n",
    "#Placeholders for the input array and the initial centroids\n",
    "input_placeholder = tf.placeholder(dtype=tf.float32, shape=[None, input_size])\n",
    "initial_centroids_placeholder = tf.placeholder(dtype=tf.float32, shape=[input_size, number_centroids])\n",
    "\n",
    "#Matrix containing the centorids\n",
    "kmeans_matrix = tf.Variable(tf.random_uniform(shape=[input_size, number_centroids], \n",
    "                                              minval=0.0, maxval=1.0, dtype=tf.float32))\n",
    "assign_centroids_op = kmeans_matrix.assign(initial_centroids_placeholder)\n",
    "\n",
    "#Here the distance is estimated and the BMU computed\n",
    "difference = tf.expand_dims(input_placeholder, axis=1) - tf.expand_dims(tf.transpose(kmeans_matrix), axis=0)\n",
    "euclidean_distance = tf.norm(difference, ord='euclidean',axis=2) #shape=(?, 3)\n",
    "bmu_index = tf.argmin(euclidean_distance, axis=1) #get the index of BMU\n",
    "bmu = tf.gather(kmeans_matrix, indices=bmu_index, axis=1) #take the centroinds\n",
    "\n",
    "#To minimize: within-cluster sum of squares (reconstruction error)\n",
    "loss = tf.reduce_mean(tf.pow(input_placeholder - tf.transpose(bmu), 2))\n",
    "\n",
    "#optimizer = tf.train.GradientDescentOptimizer(learning_rate=0.01)\n",
    "optimizer = tf.train.AdamOptimizer(learning_rate=0.05)\n",
    "\n",
    "#Training operation\n",
    "train_op = optimizer.minimize(loss=loss, global_step=tf.train.get_global_step())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now the graph is ready and we can start a session to try the algorithm. Here I simply create a certain number of random input arrays. The initialization of the centroids is trivial, I take some of the input arrays and I pass them to the centroids."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "sess = tf.Session()\n",
    "sess.run(tf.global_variables_initializer())\n",
    "input_array = np.random.random((batch_size,input_size))\n",
    "print(\"\\ninput_array.T\")\n",
    "print(input_array.T)\n",
    "print(\"\\nkmeans_matrix before assignment\")\n",
    "print(sess.run([kmeans_matrix]))\n",
    "sess.run([assign_centroids_op], \n",
    "         {initial_centroids_placeholder: input_array[0:number_centroids,:].T})\n",
    "print(\"\\nkmeans_matrix after assignment\")\n",
    "print(sess.run([kmeans_matrix]))\n",
    "\n",
    "print(\"\\nStarting training...\")\n",
    "print_every = 100\n",
    "for i in range(1000):\n",
    "    output = sess.run([train_op, loss], {input_placeholder: input_array})\n",
    "    if(i%print_every==0):\n",
    "        print(\"Loss: \" + str(output[1]))\n",
    "\n",
    "print(\"\\nkmeans_matrix final\")\n",
    "print(sess.run([kmeans_matrix]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Training on the Iris dataset\n",
    "--------------------------------\n",
    "\n",
    "Now we can try the algorithm on a real dataset. We can use the [Iris Flower dataset](../iris/iris.ipynb) for instance. The dataset is represented by the dimensions of the flowers: sepal length, sepal width, petal length, petal width. There are a total of three classes (0=Setosa, 1=Versicolor, 2=Virginica).\n",
    "The TFRecord files for this dataset are already included in Tensorbag and we can load them in memory with the following snippet:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import tensorflow as tf\n",
    "import numpy as np\n",
    "\n",
    "np_dataset = np.loadtxt(open(\"../iris/iris_dataset.csv\", \"rb\"), delimiter=\",\", skiprows=1)\n",
    "np.random.shuffle(np_dataset) #shuffle the dataset rows\n",
    "train_features_array = np_dataset[0:100, 0:4]\n",
    "train_labels_array = np_dataset[0:100, 4:].flatten().astype(np.int32)\n",
    "test_features_array = np_dataset[100:150, 0:4]\n",
    "test_labels_array = np_dataset[100:150, 4:].flatten().astype(np.int32)\n",
    "\n",
    "print(train_features_array.shape)\n",
    "print(train_labels_array.shape)\n",
    "print(test_features_array.shape)\n",
    "print(test_labels_array.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "input_size = 4 #sepal length, sepal width, petal length, petal width\n",
    "number_centroids = 3 #the number of centroids is equal to the number of the classes\n",
    "\n",
    "#Placeholders for the input array and the initial centroids\n",
    "input_placeholder = tf.placeholder(dtype=tf.float32, shape=[None, input_size])\n",
    "initial_centroids_placeholder = tf.placeholder(dtype=tf.float32, shape=[input_size, number_centroids])\n",
    "\n",
    "#Matrix containing the centorids\n",
    "kmeans_matrix = tf.Variable(tf.random_uniform(shape=[input_size, number_centroids], \n",
    "                                              minval=0.0, maxval=1.0, dtype=tf.float32))\n",
    "assign_centroids_op = kmeans_matrix.assign(initial_centroids_placeholder)\n",
    "\n",
    "#Here the distance is estimated and the BMU computed\n",
    "difference = tf.expand_dims(input_placeholder, axis=1) - tf.expand_dims(tf.transpose(kmeans_matrix), axis=0)\n",
    "euclidean_distance = tf.norm(difference, ord='euclidean',axis=2) #shape=(?, 3)\n",
    "bmu_index = tf.argmin(euclidean_distance, axis=1) #get the index of BMU\n",
    "bmu = tf.gather(kmeans_matrix, indices=bmu_index, axis=1) #take the centroids\n",
    "\n",
    "#To minimize: within-cluster sum of squares (reconstruction error)\n",
    "loss = tf.reduce_mean(tf.pow(input_placeholder - tf.transpose(bmu), 2))\n",
    "\n",
    "#optimizer = tf.train.GradientDescentOptimizer(learning_rate=0.01)\n",
    "optimizer = tf.train.AdamOptimizer(learning_rate=0.05)\n",
    "#optimizer = tf.train.RMSPropOptimizer(learning_rate=0.1)\n",
    "\n",
    "#Training operation\n",
    "train_op = optimizer.minimize(loss=loss, global_step=tf.train.get_global_step())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "tot_epochs = 1000\n",
    "\n",
    "sess = tf.Session()\n",
    "sess.run(tf.global_variables_initializer())\n",
    "print(\"\\nkmeans_matrix before assignment\")\n",
    "print(sess.run([kmeans_matrix]))\n",
    "input_array = train_features_array\n",
    "sess.run([assign_centroids_op], \n",
    "         {initial_centroids_placeholder: input_array[0:number_centroids,:].T})\n",
    "print(\"\\nkmeans_matrix after assignment\")\n",
    "print(sess.run([kmeans_matrix]))\n",
    "\n",
    "print(\"\\nLoss before training\")\n",
    "output = sess.run([loss], {input_placeholder: input_array})\n",
    "print(output)\n",
    "\n",
    "print(\"\\nStarting training...\")\n",
    "\n",
    "for iteration in range(tot_epochs):\n",
    "        output = sess.run([train_op, loss], {input_placeholder: input_array})\n",
    "        print(\"Loss: \" + str(output[1]))\n",
    "        \n",
    "print(\"\\nkmeans_matrix final\")\n",
    "print(sess.run([kmeans_matrix]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Testing the algorithm\n",
    "-------------------------\n",
    "\n",
    "After the training we can test the results. For the testing phase I assign each element to the closest centroid and I evaluate the **purity**. the purity is defined as the sum of the elements in each cluster that belong to the most numerous class, divided by the total number of elements. For instance, in the Iris dataset a purity of 1.0 can be obtained if the elements belonging to the three different classes are confined in three different clusters. It is possible to have purity one in the edge case of $K \\geq N$ meaning that the number of clusters are more than the number of elements."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#Get the indices of the test array\n",
    "indices = sess.run([bmu_index], {input_placeholder: test_features_array})[0]\n",
    "\n",
    "print(\"\\nTest classes:\")\n",
    "print(test_labels_array)\n",
    "print(\"\\nBMU indices:\")\n",
    "print indices\n",
    "\n",
    "cluster_class_matrix = np.zeros((number_centroids,3))\n",
    "for i in range(len(indices)):\n",
    "    cluster = indices[i]\n",
    "    label = int(test_labels_array[i])\n",
    "    cluster_class_matrix[cluster,label] += 1\n",
    "purity = np.sum(np.max(cluster_class_matrix, axis=1)) / np.sum(cluster_class_matrix)\n",
    "\n",
    "print(\"\\nCluster-class matrix:\")\n",
    "print cluster_class_matrix\n",
    "\n",
    "print(\"\\nPurity: \" + str(purity))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using Adam or RMSProp optimizer it is possible to obtain purity of 0.8 using three clusters. I invite you to play with the number of clusters and the optimizer to see how the purity change."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Copyright (c) 2018** Massimiliano Patacchiola, MIT License"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
